1. https://github.com/dlab-geo/r-geospatial-workshop

• Click Clone or Download and download the zip file
• Upzip the zip file and make a note of the folder in which it is located

1. Open RStudio and start a new script

2. Follow along by opening r-geospatial-workshop-pt1-sp2018.html

3. Install required libraries

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

• Intro to geospatial data

• Intro to coordinate reference systems (CRS)

• Classes & methods for working with spatial data in R

• Mapping geospatial data

• Practice

Patty Frontiera

dlab.berkeley.edu

Who are you?

Why are you here?

1. Get the workshop files

https://github.com/dlab-geo/r-geospatial-workshop

• Click Clone or Download and download the zip file

• Upzip the zip file and make a note of the folder in which it is located

1. Open RStudio and start a new script

2. Follow along by opening r-geospatial-workshop-pt1-sp2018.html

1. In RStudio, install required libraries

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

are data about locations on or near the surface of the Earth.

convey geographic information but don't unambigously specify location

represent location geometrically with coordinates

46.130479, -117.134167

Coordinates indicate specific locations on the Earth when associated with a geographic coordinate reference system or CRS.

specify

1. the shape of the Earth as the major & minor axes of an ellipsoid
2. the origin (0,0) - equator (~X), prime meridian (~Y)
3. fit of the CRS to the Earth (center of the earth), aka Datum
4. units, eg lat/lon expressed as decimal degrees or DMS

Because of variations in 1-3, there are many geographic CRSs!

The World Geodetic System of 1984 is the default geographic CRS used today.

WGS84 is the default CRS for most GIS software

Almost all data with lon/lat coordinates are assumed to be WGS84 unless otherwise specified

[NAD83] is another common geographic CRS used by US agencies like the Census.

WGS84 and NAD83 are so similar that differences are often ignored except for applications requiring high locational accuracy.

A Projected CRS applies a map projection to a Geographic CRS

Map Projection: mathematial transformation from curved to flat surface

There are many, many, many projected CRSs

All introduce distortion, eg in shape, area, distance, direction

No best one for all purposes

Selection depends on location, extent and purpose

Spatial data is a more generic term that is not just for geographic data.

• dynamically determine spatial metrics like area and length,
• characteristics like distance and direction, and
• relationships like inside and intersects from these data.

Points, lines and Polygons

Regular grid of cells (or pixels)

• We won't be covering Raster data in this workshop*

• Why special software?

Most software can only query and analyze numbers and text

• Databases, spreadsheets, statistical software packages

Software that can import, create, store, edit, visualize and analyze geospatial data

• data types/objects to represent geospatial data
  • locations represented with coordinates
  • and referenced to the surface of the earth via CRSs
• methods to operate on those representations

Desktop GIS - ArcGIS, QGIS

Software extended to support geospatial data - Tableau, R

Spatial Databases - Postgresql

Web GIS - ArcGIS Online, CARTO

Custom software - leaflet web maps

You already use R

Reproducibility

Free & Open Source

Cutting edge

Thousands of Cool add-ons

• Shiny, rleaflet

Vector points, lines & polygons:

• CSV
• ESRI Shapefile
• GeoJSON
• KML

Raster grids

• TIFF, GeoTiff
• JPEG, JPEG2000
• ASCII Grid

See GIS file formats

This is one of the most, if not the most common spatial vector data file formats.

Old but everywhere!

Gotchas: 2GB limit, 10 char column names, poor unicode support

The simplest, most common format for point data

"ID","name”,”x”,”y”,”taste","price","crowded","food"
"1","babette",-122.255374,37.868428,10,4,1,1
"2","musical",-122.260698,37.868383,7,3.25,1,1
"3","starbucks",-122.266057,37.870441,6,2.95,1,0
"4","yalis",-122.266385,37.873528,7,2.95,0,0
"5","berkeleyesp",-122.268681,37.873664,3,3.25,1,0
"6","fertile",-122.268863,37.874934,5,3.25,0,1

There are many approaches to and packages for working with geospatial data in R.

One approach is to keep it simple and store geospatial data in a data frame.

This approach is most common when the data are point data in CSV files.

San Francisco Open Data Portal https://data.sfgov.org

SF Property Tax Rolls

This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.

We are using this as a proxy for home values.

## Warning: package 'dplyr' was built under R version 3.4.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

sfhomes <- read.csv('data/sf_properties.csv')  
head(sfhomes,6)

class(sfhomes) # what type of data object?
dim(sfhomes) # how many rows and columns
str(sfhomes) # display the structure of the object
head(sfhomes) # take a look at the first 10 records
summary(sfhomes) # explore the range of values
hist(sfhomes$totvalue)  # plot he range of values for the totvalue column

plot(sfhomes$lon, sfhomes$lat) # using base plot function

library(ggplot2)

ggplot() + geom_point(data=sfhomes, aes(lon,lat), col="red", size=1)

ggplot() + geom_point(data=sfhomes, aes(lon,lat, col=totvalue))

for SalesYear equal to 2015

sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

nrow(sfhomes15) # How many records?

## [1] 889

hist(sfhomes15$totvalue) # What is the distribution of totvalue?

of sfhomes15 using ggplot and set the color to the totvalue column.

ggplot() + geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

Create basemaps on which you can display your data.

Geocode place names and addresses to get point coordinates.

and more…

library(ggmap)

# fetch map data (default=Google) to plot
# ?get_map
sf_map <- get_map("San Francisco, CA")  

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=San+Francisco,+CA&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=San%20Francisco,%20CA&sensor=false

ggmap(sf_map)

ggmap(sf_map) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

Create a basemap zoomed to the extent of our data

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes15$lon), lat = mean(sfhomes15$lat))
sf_ctr  # take a look

##        lon        lat 
## -122.43209   37.75946

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=37.759457,-122.432087&zoom=12&size=640x640&scale=1&maptype=terrain&language=en-EN&sensor=false

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

See ?get_map to get available options

sf_basemap_lite <- get_map(sf_ctr, zoom=12, scale=1, 
                            maptype = "toner-lite", source="stamen")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=37.759457,-122.432087&zoom=12&size=640x640&scale=2&maptype=terrain&sensor=false

## Map from URL : http://tile.stamen.com/toner-lite/12/653/1582.png

## Map from URL : http://tile.stamen.com/toner-lite/12/654/1582.png

## Map from URL : http://tile.stamen.com/toner-lite/12/655/1582.png

## Map from URL : http://tile.stamen.com/toner-lite/12/656/1582.png

## Map from URL : http://tile.stamen.com/toner-lite/12/653/1583.png

## Map from URL : http://tile.stamen.com/toner-lite/12/654/1583.png

## Map from URL : http://tile.stamen.com/toner-lite/12/655/1583.png

## Map from URL : http://tile.stamen.com/toner-lite/12/656/1583.png

## Map from URL : http://tile.stamen.com/toner-lite/12/653/1584.png

## Map from URL : http://tile.stamen.com/toner-lite/12/654/1584.png

## Map from URL : http://tile.stamen.com/toner-lite/12/655/1584.png

## Map from URL : http://tile.stamen.com/toner-lite/12/656/1584.png

ggmap(sf_basemap_lite) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

# Let's look at last 5 years
sfhomes2010_2015 <- subset(sfhomes, as.numeric(SalesYear) > 2009)

ggmap(sf_basemap_lite) +
  geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 )  +
  facet_wrap(~ SalesYear)

ggmap(sf_basemap_lite) +
  geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 )  +
  facet_wrap(~ SalesYear)

Redo above facet map with the following changes:

• Use data from 1995 - 1999
• Use a different basemap (eg maptype="terrain")

There are limits to what you can do with geospatial data stored in a dataframe

and mapping the data with ggplot and ggmap

Shapefile is the most common format for non-point geospatial data.

Can't directly read or plot a shapefile with data frames/ggplot/ggmap

• What properties are within walking distance (.25 miles) of a BART station?

• Average sales price by census tract?

Most spatial operations require spatial objects and methods

Classes and Methods for Spatial Data

The SP package is most commonly used to construct and manipulate spatial data objects in R.

Hundreds of other R packages that do things with spatial data typically build on SP objects.

The sf, or simple features package in R has many improvements

Based on open standards for specifying spatial data

ggplot and ggmap can map sf spatial objects

But most spatial packages still depend on sp

So, live on the bleeding edge or check back in a year or so.

Take a look at the different types of spatial object classes supported by sp

library(sp)
getClass("Spatial") 

## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

```

We use the S*DF objects most frequently!

Let's transform the sfhomes15 data frame to an sp object of type SpatialPointsDataFrame

Use the sp::coordinates() method - Sets or retrieves spatial coordinates

When transforming a DF to SPDF, requires - the object that will get the coordinates - the names of the columns that contain the X and Y coordinates

# First make a copy of the data frame
sfhomes15_sp <- sfhomes15

coordinates(sfhomes15_sp) <- c('lon','lat') # ORDER MATTERS!!

class(sfhomes15_sp) # check it

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

You transformed a data frame to an SPDF using:

coordinates(sfhomes15_sp) <- c('lon','lat')

Now try this:

coordinates(sfhomes15_sp)

str(sfhomes15) # the data frame

## 'data.frame':    889 obs. of  19 variables:
##  $ FiscalYear      : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ SalesDate       : Factor w/ 6639 levels "1920-01-30","1920-01-31",..: 6412 6507 6508 6599 6422 6536 6398 6406 6586 6421 ...
##  $ Address         : Factor w/ 166180 levels "#200 0530-SANCHEZ             ST0000",..: 40288 30412 18494 80285 143414 89595 134736 147862 2482 80378 ...
##  $ YearBuilt       : int  1978 1950 2004 1951 2015 2008 1983 1951 1963 1925 ...
##  $ NumBedrooms     : int  3 2 1 0 2 1 1 0 1 2 ...
##  $ NumBathrooms    : int  3 2 1 0 2 1 1 0 1 2 ...
##  $ NumRooms        : int  7 5 0 4 5 3 3 0 3 5 ...
##  $ NumStories      : int  2 1 0 1 1 0 1 0 0 3 ...
##  $ NumUnits        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ AreaSquareFeet  : int  2739 1147 740 836 1272 745 753 1434 775 1264 ...
##  $ ImprovementValue: int  411176 396556 446202 135000 817276 410000 85235 187312 412500 483005 ...
##  $ LandValue       : int  959410 925298 446202 315000 817276 410000 127852 561940 412500 483005 ...
##  $ Neighborhood    : Factor w/ 41 levels "","Bayview Hunters Point",..: 13 36 7 29 20 35 20 25 5 27 ...
##  $ Location        : Factor w/ 124331 levels "(37.7081643340418, -122.451068309854)",..: 73331 28165 106031 13630 90593 102763 88072 6855 115848 93189 ...
##  $ SupeDistrict    : int  7 7 6 9 8 6 8 11 3 1 ...
##  $ totvalue        : int  1370586 1321854 892404 450000 1634552 820000 213087 749252 825000 966010 ...
##  $ SalesYear       : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ lat             : num  37.8 37.7 37.8 37.7 37.8 ...
##  $ lon             : num  -122 -122 -122 -122 -122 ...

str(sfhomes15_sp) # the SPDF

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  889 obs. of  17 variables:
##   .. ..$ FiscalYear      : int [1:889] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   .. ..$ SalesDate       : Factor w/ 6639 levels "1920-01-30","1920-01-31",..: 6412 6507 6508 6599 6422 6536 6398 6406 6586 6421 ...
##   .. ..$ Address         : Factor w/ 166180 levels "#200 0530-SANCHEZ             ST0000",..: 40288 30412 18494 80285 143414 89595 134736 147862 2482 80378 ...
##   .. ..$ YearBuilt       : int [1:889] 1978 1950 2004 1951 2015 2008 1983 1951 1963 1925 ...
##   .. ..$ NumBedrooms     : int [1:889] 3 2 1 0 2 1 1 0 1 2 ...
##   .. ..$ NumBathrooms    : int [1:889] 3 2 1 0 2 1 1 0 1 2 ...
##   .. ..$ NumRooms        : int [1:889] 7 5 0 4 5 3 3 0 3 5 ...
##   .. ..$ NumStories      : int [1:889] 2 1 0 1 1 0 1 0 0 3 ...
##   .. ..$ NumUnits        : int [1:889] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ AreaSquareFeet  : int [1:889] 2739 1147 740 836 1272 745 753 1434 775 1264 ...
##   .. ..$ ImprovementValue: int [1:889] 411176 396556 446202 135000 817276 410000 85235 187312 412500 483005 ...
##   .. ..$ LandValue       : int [1:889] 959410 925298 446202 315000 817276 410000 127852 561940 412500 483005 ...
##   .. ..$ Neighborhood    : Factor w/ 41 levels "","Bayview Hunters Point",..: 13 36 7 29 20 35 20 25 5 27 ...
##   .. ..$ Location        : Factor w/ 124331 levels "(37.7081643340418, -122.451068309854)",..: 73331 28165 106031 13630 90593 102763 88072 6855 115848 93189 ...
##   .. ..$ SupeDistrict    : int [1:889] 7 7 6 9 8 6 8 11 3 1 ...
##   .. ..$ totvalue        : int [1:889] 1370586 1321854 892404 450000 1634552 820000 213087 749252 825000 966010 ...
##   .. ..$ SalesYear       : int [1:889] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   ..@ coords.nrs : int [1:2] 19 18
##   ..@ coords     : num [1:889, 1:2] -122 -122 -122 -122 -122 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:889] "151677" "150126" "125377" "65555" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] -122.5 37.7 -122.4 37.8
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

You can see from str(sfhomes) that a SPDF object is a collection of slots or components. The key ones are:

• @data data frame of attributes that describe each location
• @coords the coordinates for each geometric object - here points
• @bbox the min and max bounding coordinates
• @proj4string the coordinate reference system defintion as a string

Explore the object in the Environment window

Review the output of each of these:

summary(sfhomes15_sp)
head(sfhomes15_sp@data)
class(sfhomes15_sp@data)

sfhomes15_sp@bbox
bbox(sfhomes15_sp)

head(sfhomes15_sp@coords)
head(sfhomes15_sp$lat)
head(sfhomes15_sp$lon)

sfhomes15_sp@proj4string
proj4strings(sfhomes15_sp)

Look at sfhomes15_sp in the environment window

Are all the columns that were present in sfhomes15 also in sfhomes15_sp?

Is there a slot in sfhomes15_sp without data?

proj4string(sfhomes15_sp) # get a CRS object

## [1] NA

plot(sfhomes15_sp)  # using sp::plot

We created great maps of sfhomes point data with ggplot and ggmap.

Then we created a simple map of the SPDF with plot.

We aren't seeing the value of sp objects just yet.

There's an R package for that!

rgdal is an R port of the powerful and widely used GDAL library.

It is the most commonly used R library for importing and exporting spatial data.

• OGR: for vector data: readOGR() and writeOGR()

• GDAL for raster data: readGDAL() and writeGDAL()

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.4.3

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-5

# See what file types are supported by rgdal drivers
# ogrDrivers()$name

gdal.org

`?readOGR

For more info on working with rgdal to load different types of spatial data in R see this excellent tutorial by Zev Ross.

Take a look at the file(s)

dir("data", pattern="sf_boundary")

## [1] "sf_boundary.dbf" "sf_boundary.prj" "sf_boundary.shp" "sf_boundary.shx"

sfboundary <- readOGR(dsn="data",layer="sf_boundary")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_boundary"
## with 1 features
## It has 3 fields

# or
# sfboundary <- readOGR("data","sf_boundary")
# but not
#sfboundary <- readOGR(dsn="data/",layer="sf_boundary")

What type of sp object is sfboundary? How many features are in the object?

class(sfboundary)
str(sfboundary) 
head(sfboundary@data)  

Explore the object in the Envi window

How?

plot(sfboundary)

library(tigris)
calcounties <- counties(state="California", cb=T)
class(calcounties)
sf_cen <- subset(calcounties, COUNTYFP == "075")
plot(sf_cen)
head(sf_cen@data)
# Now look in Envi Window

Is it a geographic or projected CRS?

proj4string(sfboundary)

plot(sfboundary)
points(sfhomes15_sp, col="red")

Where are the points? What's wrong?

plot(sfboundary)
points(sfhomes15_sp, col="red")

Compare the CRSs, are they the same?

proj4string(sfboundary)
proj4string(sfhomes15_sp)
proj4string(sfboundary) == proj4string(sfhomes15_sp)

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] NA

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] NA

sfboundary@bbox

##         min       max
## x  542696.6  556659.9
## y 4173563.7 4185088.6

sfhomes15_sp@bbox

##            min        max
## lon -122.51065 -122.36964
## lat   37.70821   37.80586

The #1 reason…

All sp objects should have a defined CRS

If not, one needs to be assigned to the object - This is called defining a projection. - This doesn't change the coordinates.

All sp objects should have the same CRS.

• When they don't, they need to be transformed to a common CRS.

• This is also called a projection transformation,
  • or projecting or reprojection.

• Projection transformation returns a new spatial object with the transformed coordinates

Geospatial-aware software will have a database of definitions for thousands of Earth referenced coordinate systems

We need to define the CRS of what sp object? - sfboundary or sfhomes15_sp?

We need to know - the appropriate CRS for the data - how to define the CRS - how to assign it to the sp object

sp includes the CRS() function to define a CRS

• Two ways

# use an EPSG code
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

# or enter the proj4 string
# proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
#                               +ellps=WGS84 +datum=WGS84 +no_defs")  

proj4string(sfhomes15_sp)

## [1] "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

• proj4string() can get or set the CRS
• proj4 is a library for managing map projections and CRSs
• epsg codes are used as short-hand for CRSs definitions

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] FALSE

Use sp function spTransform

Requires as input:

• a sp object to transform with a defined CRS

• a target CRS

Outputs a new spatial object with coordinate data in the target CRS

• It does not change the input data

from UTM10 CRS to WGS84 (why?)

sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))

# or
# sfboundary_lonlat <- spTransform(sfboundary, 
#                              CRS(proj4string(sfhomes15_sp)))

How do these two spTransform approaches differ?

How will we know?

proj4string(sfhomes15_sp) == proj4string(sfboundary_lonlat)

plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
points(sfhomes15_sp[sfhomes15_sp$totvalue<1000000,], col="green")

Woo-hoo!

We can transform sp objects to the same CRS for mapping and spatial analysis (we'll get there!)

Use writeOGR to save sfboundary_lonlat to a new shapefile

See ?writeOGR for help

# write transformed data to a new shapefile 
writeOGR(sfboundary_lonlat, 
          dsn = "data", 
          layer = "sfbounary_lonlat", 
          driver="ESRI Shapefile")

# is it there?
dir("data")

We want all data in the same CRS

Which one is best?

See http://spatialreference.org/

Use this site to find EPSG codes and proj4 CRS strings

Geographic CRSs

• 4326 Geographic, WGS84 (default for lon/lat)

• 4269 Geographic, NAD83 (USA Fed agencies like Census)

Projected CRSs

• 5070 USA Contiguous Albers Equal Area Conic

• 3310 CA ALbers Equal Area

• 26910 UTM Zone 10, NAD83 (Northern Cal)

• 3857 Web Mercator (web maps)

base::plot, ggplot, ggmap for geospatial data in data frames

• great for creating maps given these types of data

sp::plot for sp objects

• meh, but great for a quick look at spatial data

Use it to create quick maps

Can be used to create great maps

• BUT complex, non-intuitive syntax = long ugly code

See examples in sp Gallery: Plotting maps with sp

# map of the sfhomes data by totalvaue
spplot(sfhomes15_sp,"totvalue")

tmap stands for thematic map

Great maps with less code than the alternatives

Syntax should be familar to ggplot2 users, but simpler

Relatively easy to create interactive maps

• ?tmap
• vignette("tmap-nutshell")
• vignette("tmap-modes")
• tmap in a Nutshell

Load the library

library(tmap)

## Warning: package 'tmap' was built under R version 3.4.3

qtm(sfhomes15_sp)

tmap_mode("view")

## tmap mode set to interactive viewing

qtm(sfhomes15_sp)

tmap_mode("plot")

## tmap mode set to plotting

NOTE: columnn names must be quoted!

qtm(sfhomes15_sp, symbol.col="totvalue")

Create a static qtm of sfboundary_lonlat

• bonus, set the border color to black and the fill to beige

See ?qtm for description of all options

qtm(sfboundary_lonlat, borders="black", fill="beige")

tmap's flexibility comes in how it intitively allows you to layer spatial data and style the layers by data attributes

Use tm_shape(<sp_object>) to specifiy a geospatial data layer

Add + tm_<element>(...) to style the layer by data values

…and other options for creating a publication ready map

?tmap_shape

?tmap_element

• ?tm_polygons (tm_fill, tm_borders)

• ?tm_symbols (tm_dots, etc…)

• tmap in a Nutshell

tm_shape(sfboundary_lonlat) + tm_polygons(col="beige", border.col="black")

tm_shape(sfhomes15_sp) + tm_dots(col="totvalue", size=.5)

Make the previous map interactive

Zoom in on light & dark colored points to see if the corressponding totvalue makes sense.

tmap_mode("plot")

## tmap mode set to plotting

tm_shape(sfboundary_lonlat) + tm_polygons(col="black", border.col="grey") + 
   tm_shape(sfhomes15_sp) + tm_dots(col="totvalue", size=.5)

Redo the last map but

change the legend title to "San Francisco Home Sales, 2015"

• hint: title=

We have been mapping data in WGS84 (lon/lat) CRS.

Not a great idea for static maps of larger areas as distortion becomes evident.

Let's explore that distortion and how to address with tmap

Let's load and map data for US states.

Data are in the file data/us_states_pop.shp

Use readOGR to load data/us_states_pop.shp into an sp object named us_states

us_states <- readOGR("./data", "us_states_pop")

## OGR data source with driver: ESRI Shapefile 
## Source: "./data", layer: "us_states_pop"
## with 49 features
## It has 4 fields

qtm(us_states)

Review us_states with the sp commands we used earlier and / or explore in the Environment window.

• What type of sp object is us_states

• How many features does it contain?

• How many attributes describe those features?

• What is the CRS?

tm_shape(us_states) + tm_polygons(col="grey", border.col = "white")

Notice anything odd about shape of USA?

tm_shape(us_states, projection="+init=epsg:5070") + tm_polygons(col="grey", border.col = "white")

tm_shape(us_states, projection="+init=epsg:5070") + tm_polygons(col="grey", border.col = "white") +
tm_shape(us_states) + tm_borders(col="purple") 

Also called On-the-fly reprojection in ArcGIS & QGIS

Very cool!

BUT, if you want to use data in a different CRS it is best to transform it.

• Transform the us_states data to the USA Contiguous Albers CRS (5070),

• Save output as a new SpatialPolygonsDataFrame called us_states_5070

us_states_5070 <- spTransform(us_states, CRS("+init=epsg:5070"))

tm_shape(us_states_5070) + tm_polygons(col="beige") +
  tm_shape(us_states) + tm_borders(col="purple")